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ABSTRACT 

We describe S2LET, a fast and robust implementation of the scale-discretised wavelet transform on the sphere. Wavelets are con- 
structed through a tiling of the harmonic line and can be used to probe spatially localised, scale-depended features of signals on the 
sphere. The scale-discretised wavelet transform was developed previously and reduces to the needlet transform in the axisymmetric 
case. The reconstruction of a signal from its wavelets coefficients is made exact here through the use of a sampling theorem on the 
sphere. Moreover, a multiresolution algorithm is presented to capture all information of each wavelet scale in the minimal number of 
samples on the sphere. In addition S2LET supports the HEALPix pixelisation scheme, in which case the transform is not exact but 
nevertheless achieves good numerical accuracy. The core routines of S2LET are written in C and have interfaces in Mat lab, IDL 
and Java. Real signals can be written to and read from FITS files and plotted as Mollweide projections. The S2LET code is made 
publicly available, is extensively documented, and ships with several examples in the four languages supported. At present the code 
is restricted to axisymmetric wavelets but will be extended to directional, steerable wavelets in a future release. 
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1. Introduction 

Signals defined or measured on the sphere arise in numerous 
disciplines, where analysis techniques defined explicitly on the 
sphere are now in common use. In particular, wavelets on the 
sphere (|Antoine & V andergh eynst|1998][T999l I Baldi et al.|2009 
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2006||Starck et al.|2(X)6a||Wiaux et al.|2006||2005||2007[ T2008 1 
Yeo et al. 2008 ) have been applied very successfully to problems 
in astrophysics and cosmology, where data-sets are increasingly 
large and need to be analysed at high resolution in order to con- 
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While wavelet theory is well established in Euclidean space 
(see e.g. Daubechies (1992)), multiple wavelet frameworks have 
been developed on the sphere, only a fraction of which lead to 
exact transforms in both the continuous and discrete settings. 
In fact, discrete methodologies ( |Schroder & Sweldens| [l995; 
Sweldens 1996 1997| ) achieve exactness in practice but may not 



lead to a stable basis on the sphere ( Sweldens|1997| l. In the con- 
tinuous setting several constructions are theoretically exact, and 
have been combined with sampling theorems on the sphere to 
enable exact reconstruction in the discrete setting also. In par- 
ticular, scale-discretised wavelets (Wi aux et al.|20 08) lean on a 



tiling of the harmonic line to yield an exact wavelet transform 
in both the continuous and discrete settings. In the axisymmetric 
case, the scale-discretised wavelets (Wiaux et al. 2008) reduce 
to needlets ( |Narcowich et al.|2006||Baldi et al.|2009[|Marinucci| 
et al. 2008} , which were developed independently using a similar 
tilling of the harmonic line. 

In this paper we describe the new publicly available S2LElQ 
code to perform the scale-discretised wavelet transform of com- 
plex signals on the sphere. At present S2LET is restricted to 
axisymmetric wavelets (i.e. azimuthally symmetric when cen- 
tred on the poles) and includes generating functions for both ax- 



isymmetric scale-discretised wavelets (Wiaux et al. 
needlets (jNarcowich et al.||2006l [Baldi et al.||2009| 



[et~aT1 [2008 



2008 ) and 



Marinucci 



We intend to extend the code to directional, 
steerable wavelets and spin functions in a future release. The 
core routines of S2LET are written in C, exploit fast algo- 
rithms on the sphere, and have interfaces in Mat lab, IDL and 
Java. We note that public cod es are already ava ilable to com- 
pute scale-discr etised (S2DIaF1 [Wiaux et aL 2008 1 and needlet 
(NeedATooIp| |Pietrobon et al'|2010| ) transforms. We envisage 
S2LET will supersede S2DW once it is extended to directional, 
steerable wavelets and spin functions. 

The remainder of this article is organised as follows. In sec- 
tion 2 we detail the construction of scale-discretised axisym- 
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metric wavelets and the corresponding exact scale-discretised 
wavelet transform on the sphere. In section 3 we describe the 
S2LET code, including implementation details, computational 
complexity and numerical performance. We present a number of 
simple examples using S 2 LET in section 4, along with the code 
to execute them. We conclude in section 5. 



2. Wavelets on the sphere 

We review scale-discretised wavelets on the sphere, as consid- 
ered previously by Wi aux et al.| (2008i. Directional, steerable 
wavelets were also considered by Wiaux et al.| ([2008 ), how- 
ever we restrict our attention to axisymmetric wavelets here. The 
use of a sampling theorem on the sphere guarantees that spher- 
ical harmonic coefficients capture all the information content of 
band-limited signals, resulting in exact harmonic and wavelet 
transforms; hence, we recall harmonic analysis on the sphere 
concisely. One may alternatively adopt samplings of the sphere 
for which exact quadrature rules do not exist such as HEALPix 
( |G6rski et al.|2005] >, but which nevertheless exhibit other useful 
properties, leading to numerically accurate but not theoretically 
exact transforms. 



2.1. Harmonic analysis on the sphere 

The spherical harmonic decomposition of a square integrable 
signal / e L 2 (S 2 ) on the two-dimensional sphere S 2 reads 



(1) 



e=o m=-l 



where Y( m are the spherical harmonic functions, which form the 
canonical orthogonal basis on S 2 . The spherical harmonic coef- 
ficients f( m , with t G N such that I < oo and m e Z such that 
\m\ < t, form a dual representation of the signal / in the har- 
monic basis on the sphere. The angular position u> — (0, <p) e S 2 
is specified by colatitude 6 e [0, n] and longitude (f> e [0, 2n). 
The spherical harmonic coefficients are given by 



f(m - if\Y{ m ) 



-L 



dn( w )/( w )y;,M, 



(2) 



with the surface element dQ.(oj) = sm6d9d<p. The signal / is 
band-limited in the spherical harmonic basis with band-limit L 
if ftm — 0, W > L. For band-limited signals sampling theo- 
rems can be invoked so that both forward and inverse transforms 
can be reduced to finite summations that are theoretically exact. 
Sampling theorems effectively encode a quadrature rule for the 
exact evaluation of integrals on the sphere from a finite set of 
sampling nodes. Various sampling theorems exist in the litera- 
ture (e.g. I Driscoll & He aly 1994; |Healy et al.|[T99"6} |McEwen 
|& Wiaux] 201 l| l. In this work we adopt the |McEwen & W iaux 
( 201 l| l sampling theorem (hereafter MW), which is based on an 
equiangular sampling scheme and, for a given band-limit L, re- 
quires the lowest number of samples on the sphere of all sam- 
pling theorems, namely (L - 1)(2L - 1) + 1 ~ 2L 2 samples (for 
comparison ~ 4L 2 samples are required by |Driscoll & Healy 
( |1994| l). Fast algorithms to compute the corresponding spherical 
harmonic transform scale as 0(L 3 ) and are numerically stable 
to band-limits of at least L = 4096 ( |McEwen & Wiaux ||20 1 1) . 
Alternative sampling schemes that are not based on sampling 
theorems also exist, such as HEALPix ( |G6rski et al.||2005] l or 
GLESP ( |Doroshkevich et aT|2 005). However, these schemes do 
not lead to exact spherical harmonic transforms on the sphere. 



Nevertheless the resulting approximate transforms achieve good 
accuracy and benefit from other practical advantages, such as 
equal-area pixels. 

2.2. Scale-discretised wavelets on the sphere 

The scale-discretised wavelet transform aims to exact spatially 
localised, scale-dependent features in the signal of interest 
/ e L 2 (S 2 ). The j-th wavelet coefficient W vJ € L 2 (S 2 ) is defined 
as the convolution of / with the wavelet e L 2 (S 2 ): 



f(6))s(f*^) = ^ 



dD.(oj')f(co')(K 1 Wr(cj'), 



(3) 



where * denotes complex conjugation. Convolution on the sphere 
is defined by the inner product of / with the rotated wavelet 
< R ll) v if-'. We restrict our attention to axisymmetric wavelets, i.e. 
wavelets that are azimuthally symmetric when centred on the 
poles. Consequently, the rotation operator is parameterised 
by angular position co = (0, (/>) only and not also orientationj^For 
the axisymmetric case the spherical harmonic decomposition of 
W^' is then simply given by a weighted product in harmonic 
space: 



ry (m 



An 
27+ 1 



. ,,„ - {W*'\Y tm ), f im = {f\Y tm ) and ^«5 m0 



(4) 



where W%, 

and where 6 m o is the Kronecker delta symbol. 

The wavelet coefficients extract the detail information of 
the signal only; a scaling function and corresponding scaling 
coefficients must be introduced to represent the low-frequency 
(low-f), approximate information of the signal. The scaling co- 
efficients W° e L 2 (S 2 ) are defined by the convolution of / with 
the scaling function <E> e L 2 (S 2 ): 



lf(w)E(f+0)(w) = (/W>, 



or in harmonic space, 



W? = 



4n 
2(+\ 



fem®*e , 



(5) 



(6) 



where Wf m = <W*|Fft,> and <t> eo 6 m0 = (<b\Ytm). 

Provided the wavelets and scaling function satisfy an admis- 
sibility property (defined below), the function / may be recon- 
structed exactly from its wavelet and scaling coefficients by 

f(oS) = ^ dQ( W ')W°(a/)(!^<<I>)(^) 



j=Jo JS 



dnrvw 'VXK/^Xw), (7) 



or equivalently in harmonic space by 



j=Jo 



4 As already noted, the extension to directional scale-discretised 
wavelets has been derived by |Wiaux et al.| p008 ). At present the S 2 LE T 
code supports axisymmetric wavelets only; directional wavelets will be 
added in a later release. 
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The parameters Jo, J define the lowest and highest scales j of 
the wavelet decomposition and must be defined consistently to 
extract and reconstruct all the information content of /. These 
parameters depend on the construction of the wavelets and scal- 
ing function and are defined explicitly in the next paragraphs. 
The admissibility condition under which a band-limited func- 
tion / can be decomposed and reconstructed exactly is given by 
the following resolution of the identity: 



An 

21 +\ 
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j=Jo 



= 1, W. 



(9) 



We are now in position to define wavelets and a scaling func- 
tion that satisfy the admissibility property. We use the smooth 
generating functions defined by Wiau x et al.| ([2008 ) in order to 
tile the harmonic line (the needlet generating functions defined 
by |Marinucci et al.| ( |2008) can also be used; in fact, the S 2 LET 
code supports both definitions). Consider the C°° Schwartz func- 
tion with compact support on [-1, 1]: 



s(t) = 



e 
0, 



t € [-1,1] 



for t e R. We introduce the positive real parameter A e 
map s(t) to 

/ 2A 



(10) 



to 



(11) 



which has compact support in [1/A, 1]. We then define the 
smoothly decreasing function k A by 



hit) = 



J 1 fftn 



(12) 



which is unity for t < I /A, zero for t > 1, and is smoothly 
decreasing from unity to zero for t e [1//1, 1], We finally define 
the wavelet generating function by 




8 , 16 32 
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(a) Tiling of the harmonic line 



64 126 




(b) Wavelets and scaling function on the sphere 

Fig. 1: Wavelets for scales j e {2,3,4,5} and scaling function 
constructed through a tiling of the harmonic line (top panel) and 
reconstructed on the sphere (bottom panel) for the parameters 
A — 3 and Jo — 2 and for band-limit L = 128. This plot was 
produced with a Mat lab demo included in S2LET. 



K A {t) = y/k A {t/A) - k A (t) 
and the scaling function generating function by 

ift(0 = Vw 



(13) 



(14) 



The wavelets and scaling function are constructed from their 
generating functions to satisfy the admissibility condition given 
by Eqn. i9i. A natural approach is to define from the gener- 
ating functions k a to have support on [ / V _1 ,/l- /+1 ], yielding 



arbitrary, provided that < Jo < J. The wavelets and scal- 
ing function may then be reconstructed on the sphere through 
an inverse spherical harmonic transform. With this construction 
the wavelets and scaling function are well-localised both spa- 
tially on the sphere and also in harmonic space, as shown in 
Figure [T] Consequently, the scale-discretised wavelet transform 
on the sphere can be used to extract spatially localised, scale- 
dependent features in the signal of interest. 



tin 
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(15) 



For these wavelets, Eqn. ^ is satisfied for t > A J °, where 7 is 
the lowest wavelet scale used in the decomposition. The scaling 
function <1> is constructed to extract the modes that cannot be 
probed by the wavelets (i.e. modes with t < A J °): 



2£+\ It, 



(16) 



To satisfy exact reconstruction, J is set to ensure the 
wavelets reach the band-limit of the signal of interest, yielding 
/ = riog,(L- 1)]. The choice of the lowest wavelet scale Jo is 



3. The S2LET code 

In this section we describe the S2LET code. We first introduce 
a multiresolution algorithm to capture each wavelet scale in the 
minimum number of samples on the sphere, which follows by 
taking advantage of the reduced band-limit of the wavelets for 
scales j < J - 1 . This multiresolution algorithm was introduced 
by Wiaux et al. ( 2008 1 already and reduces the computation cost 



of the transform considerably. We then provide details of the im- 
plementation, the computational complexity and the numerical 
accuracy of the scale-discretised wavelet transform implemented 
in the S2LET code. We finally outline planed future extensions 
of the S2LET code. 
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3.1. Multiresolution algorithm 

In harmonic space, the wavelet coefficients are simply given 
by the weighted product of the spherical harmonic coefficients 
of / and the wavelets, as expressed in Eqn. (|6Jl. Although the 
wavelet coefficients can be analysed at the same resolution as 
the signal / (full-resolution), by construction they have differ- 
ent band-limits for different scales j, as shown in Figure [T] 
Reconstruction can thus be performed at a lower resolution with- 
out any loss of information. This approach yields a multiresolu- 
tion algorithm where the wavelet coefficients are reconstructed 
with the minimal number of samples on the sphere: the y'-th 
wavelet coefficients have band-limit A J+1 and are thus recovered 
on (A j+1 - l)(2/f +1 - 1) + 1 samples on the sphere when the MW 
sampling theorem is adoptedPj This approach leads to signifi- 
cant improvements in terms of speed and memory use compared 
to the full-resolution case, as shown in the next section. Figure[2] 
illustrates the use of the full-resolution and multiresolution trans- 
forms on a map of Earth topography data. 



3.2. Implementation and interfaces 

The core numerical routines of S2LET are implemented in C. 
By adopting a low level programming language such as C for the 
implementation of the core algorithms, computational efficiency 
is optimised. The C library includes the full-resolution and 
multiresolution wavelet transforms, with specific optimisations 
for real signals in order to take advantage of all symmetries of 
the spherical harmonic transform. The scale-discretised wavelet 
transform is computed in harmonic space through Eqn. Q and 
Eqn. |6|, for the input parameters (L, A, Jo). To reconstruct sig- 
nals on the sphere, by default S2LET uses the exact spherical 
harmonic tra nsform of the MW sampling theorem (McEwen & 
Wiaux||201 1 1 implemented in the SSHlj^code. In this case all 
transforms are exact and one can analyse and synthesise real and 
complex signals using the MW sampling at floating-point preci- 
sion. S2LET has been extended to also support the HEALPix 
sampling scheme. Since the spherical harmonic transform im- 
plemented in HEALPix is not theoretically exact, the scale- 
discretised wavelet transform computed using the HEALPix 
sampling does not capture all the information of a band-limited 
signal and is thus not theoretically exact. Nevertheless, the trans- 
forms computed using HEALPix achieve good numerical accu- 
racy in most situations. 

We provide interfaces for the C library in three languages: 
Matlab, IDL and Java. The Matlab and IDL codes also 
include routines to read/write signals on the sphere stored in 
either HEALPix F Unfiles or the FITS file format used to 
stored MW sampled signals. In addition, functionality to plot the 
Mollweide projection of real signals for both MW or HEALPix 
samplings is included. The Java interface includes an object- 
oriented representation of sampled maps, spherical harmonics 
and wavelet transforms. All routines and interfaces are well doc- 
umented and illustrated with several examples for both the MW 
and HEALPix samplings. These examples cover multiple com- 
binations of parameters and types of signals. S2LET requires 
SSHT, which implements fast and exact algorithms to perform 



5 When adopting the HEALPix sampling of the sphere, the 
HEALPix resolution parameter N J side for the wavelet coefficients of 
scale j can be related to the band-limit of the wavelet coefficients: e.g. 
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(a) Full-resolution wavelet transform 




l \ / 
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(b) Multiresolution wavelet transform 

Fig. 2: Scale-discretised wavelet transform of a band-limited to- 
pography map of the Earth for A = 3, Jq = 2 and L = 128, i.e. 
with the wavelets shown in Figure [T] The wavelet transform de- 
composes the band-limited signal into wavelet coefficients that 
extract spatially localised, scale-dependent features. Since the 
wavelets for different scales j have different band-limits, the 
wavelet coefficients can be reconstructed at lower resolution on 
the sphere for lower scales /'. Panel (a) shows the full-resolution 
wavelet transform of the topography map. The original Earth to- 
pography map is shown in the top-left plot, the scaling coef- 
ficients are shown in the top-right plot, while the wavelet co- 
efficients at scales j e {2,3,4,5} are shown left-to-right, top- 
to-bottom respectively in the remaining plots. Panel (b) shows 
the same decomposition but using the multiresolution algorithm. 
The signals shown in panel (b) contain the same information as 
in panel (a) but represented in the minimal number of samples 
on the sphere. These plots were produced by one of the many 
Matlab demos provided with S2LET. 



the forward and inverse spherical harmonic transforms corre- 
sponding to the MW sampling theorem (McEwen & Wiaux 
|201l| . SSHT in turn requires the FFTt^package for the com- 
putation of fast Fourier transforms. The fast spherical harmonic 
transforms implemented in SSHT compute Wigner functions, 
and thus the spherical harmonic functions, through efficient re- 
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cursion using either the method of Trapani & Navaza (2006) 
Risbo ( 1996). Here we present results using the recursion of 



Risbo| ( |1996[ ). The fast spherical harmo nic transform algorithms 
implemented in SSHT scale as 0(L 3 ) (McE wen & Wiaux|20lT] l. 

The complexity of the axisymmetric scale-discretised 
wavelet transform is dominated by spherical harmonic trans- 
forms since the wavelet transforms are computed efficiently in 
harmonic space, through Eqn. (|4i and Eqn. (|6]l for the for- 
ward transform and through Eqn. ( 8} for the inverse transform. 
Given a band-limit L and wavelet parameters (A, Jo), recall that 
the maximum scale is given by J — \\og A (L - 1)1 and hence 
the wavelet transform (forward or inverse) involves (J - Jo + 3) 
spherical harmonic transforms (one for the original signal, one 
for the scaling coefficients and (J - Jo + 1) for the wavelet co- 
efficients). If the scaling coefficients and all wavelet coefficients 
are reconstructed at full-resolution in real space, the axisymmet- 
ric scale-discretised wavelet transform scales as 0((J-Jq+3)L?). 
However, in the previous section we established a multiresolu- 
tion algorithm that takes advantage of the reduced band-limit 
of the wavelets for scales j < J — I, With the multiresolution 
algorithm only the finest wavelet scales j e {J - 1, J} are com- 
puted at maximal resolution corresponding to the band-limit of 
the signal. The complexity of the overall multiresolution wavelet 
transform is then dominated by these operations and effectively 
scales as 0{L?). 



3.3. Numerical validation 

We evaluate the performance of S 2 LET in terms of accuracy and 
complexity using the MW sampling theorem, for which all trans- 
forms are theoretically exact. We show that S2LET achieves 
floating-point precision and scales as detailed in the previous 
section. 

We consider band-limits L — 2' with i 6 {2, . . . , 10} and gen- 
erate sets of spherical harmonic coefficients fy m following in- 
dependent Gaussian distributions 7V(0, 1). We then perform the 
wavelet decomposition and reconstruct the harmonic coeffi- 
cients, denoted by f^. We evaluate the accuracy of the trans- 
form using the error metric e = max \ff m - /fff I, which is the- 
oretically zero for both transforms since all signals are band- 
limited by construction. The complexity is quantified by observ- 
ing how the computation time t c - [f syn thesis + fanaiysisl /2 scales 
with band-limit, where the synthesis and analysis computation 
times are specified by f syn thesis and f an aiysis respectively. Since we 
evaluate the wavelet transform in real space, a preliminary step 
is required to reconstruct the signal / from the randomly gener- 
ated f[ m . This step is not included in the computation time since 
its only purpose is to generate a valid band-limited test signal 
on the sphere. The analysis then denotes the decomposition of / 
into wavelet coefficients W w ' and scaling coefficients W 4 * on the 
sphere. The synthesis refers to recovering the signal / rec from 
these coefficients. The final step, which is not included in the 
computation time either, is to decompose f ec into harmonic co- 
efficients in order to compare them with f( m . The stability of 
both e and t c is checked by averaging over hundreds of realisa- 
tions of fc, n for L-2 1 with i 6 {2, . . . , 8} and a few realisations 
with i e {9, 10). The results proved to be very stable, i.e. the vari- 
ances of the error and timing metrics are lower than 5%. Recall 
that for given band-limit L the number of samples on the sphere 
required by the exact quadrature is (2L - 1)(L - 1) + 1. All tests 
were run on an Intel 2.0GHz Core i7 processor with 8GB of 
RAM. 
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(a) Numerical accuracy of the wavelet transform 
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(b) Computation time of the wavelet transform 
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Fig. 3: Numerical accuracy and computation time of the scale- 
discretised wavelet transform computed with S2LET. We con- 
sider L — 2' with i e {2, . . . , 10), with parameters A = 2, Jo = 0. 
These results are averaged over many realisations of random 
band-limited signals and were found to be very stable. The scale- 
discretised transform is either performed at full-resolution (solid 
lines) or with the multiresolution algorithm (dashed lines). Very 
good numerical accuracy is achieved by both the full-resolution 
and multiresolution algorithms (which achieve indistinguishable 
accuracy), with numerical errors comparable to floating-point 
precision, found empirically to scale as O(L) as shown by the 
red line in panel (a). Computation time scales as <9(L 3 ) for both 
algorithms as shown by the red line in panel (b), in agreement 
with theory. The multiresolution algorithm is four to five times 
faster than the full-resolution approach for the band-limits con- 
sidered. 



The accuracy and timing performance of the scale- 
discretised wavelet transform implemented in S2LET are pre- 
sented in Figure [3] S2LET achieves very good numerical ac- 
curacy, with numerical errors comparable to floating-point pre- 
cision. Moreover, the full-resolution and multiresolution algo- 
rithms are indistinguishable in terms of accuracy. However, 
the latter is four to five times faster than the former for the 
band-limits considered since only the wavelet coefficients for 
j e {J—1,J} are computed at full-resolution. As shown in 
Figure [3] computation time scales as 0(L?) for both algorithms, 
in agreement with theory. 

3.4. Future extensions 

We have restricted this work to axisymmetric scale-discretised 
wavelets so that the rotation operator is only parametrised by 
angular position on the sphere and not also orientation. The ex- 
tension to directional wavelets has already been considered by 
|Wiaux et al.| ([2008 ) but has not yet been implemented in the 
S2LET code. In future work we plan to relax this constraint and 
to consider directional, steerable wavelets on the sphere (Wiaux 
|et al.|[2008] ). We also plan to exploit recent ideas leading to 
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fast (spin) spherical harmonic transforms (McEwen & Wiaux 



2011 



et al. 



to yield faster algorithms than those developed by Wiaux 



( 2008 ) to compute directional wavelet transforms on the 



sphere. Finally, we intend to add support to analyse spin signals 
on the sphere. In a future release, the code will also be paral- 
lelised, which will lead to further speed improvements, and a 
more complete framework in Java will be added, including IO 
and plotting routines. The S2LET code will thus be under active 
development with future releases forthcoming. In any case, we 
hope this first version of the S2LET code will prove useful for 
axisymmetric scale-discretised wavelet analysis on the sphere. 
Indeed, the code has already been used as an integral part of 
the new exact flaglet wavelet transform on the ball (Leistedt & 
McE wen|2012 1, the spherical space constructed by augmenting 
the sphere with the radial line. 



4. Examples 

The S2LET code is extensively documented and ships with sev- 
eral examples in the four languages supported. In this section we 
present a subset of short examples, along with the code to exe- 
cute them in order to demonstrate the ease of using S2LET to 
perform wavelet transforms j^] 

4. 1 . Wavelet transform from the command line 

S 2 LET includes ready-to-use high-level programs to directly de- 
compose a real signal into wavelet coefficients. The inputs are a 
FITS file containing the signal of interest and the parameters 
for the transform. The program writes the output coefficients 
in FITS files in the same directory as the input file and with 
a consistent naming scheme. These commands are available for 
both sampling schemes. For the MW sampling case illustrated in 
Example[T] the wavelet transform is exact and the band limit cor- 
responds to the resolution of the input map, which will be read 
automatically from the file. The transform may be performed in 
full-resolution or multiresolution by adjusting the multiresolu- 
tion flag specified by the last parameter (respectively and 1), 
and the output wavelet coefficients are computed at full and min- 
imal resolution accordingly. For the case of a HE ALP ix map, as 
illustrated in Example[2] the transform is not exact and the band- 
limit must be supplied as the last parameter in the command. The 
output scaling and wavelet coefficients of a HEALPix map are 
reconstructed and stored in FITS files at the same resolution as 
the input map. For both MW and HEALPix samplings the out- 
put coefficients may be read and plotted using the Mat lab or 
IDL routines. 



>> . /bin/s21et_axisym_hpx_analysis_real 
<inputFitsFile> <lambda> <J_0> 
<bandLimit> 

>> . /bin/s21et_axisym_hpx_synthesis_real 

<outputRoot> <lambda> <J_0> <bandLimit> 



Example 2: Performing the forward (analysis) and inverse 
(synthesis) wavelet transform of a real signal (HEALPix sampling) 
from the command line. 



% Example: Wavelet transform in Matlab 
lambda =3; JO = 2; L = 192; 
Jmax = s21et_jmax (L, lambda); 

% Read a real HEALPix map from a FITS file 
inputfile = ' data/somecmbsimu_hpx_12 8 . fits' ; 
[f, nside] = 

s21et_hpx_read_real_map (inputfile) ; 

% Perform the wavelet transform 
[f_wav, f_scal] = s21et_axisym_hpx_analysis 
(f , ' B' , lambda, ' L' , L, ' J_min' , JO) ; 

% Plot the map and the wavelet coefficients 
figure; ns = ceil (sqrt (2+Jmax-JO+l ) ) ; 
subplot (ns, ns, 1); 
s21et_hpx_plot_mollweide (f ) ; 
title (' Initial band-limited data') 
subplot (ns, ns, 2); 
s21et_hpx_plot_mollweide (f_scal) ; 
title (' Scaling fct') 
for j = JO : Jmax 

subplot (ns, ns, j-JO+3); 

s21et_hpx_plot_mollweide ( f_wav { j- JO+1 } ) ; 
title ([' Wavelet scale : 
' , int2str ( j) -JO+1] ) 

end 



>> . /bin/s21et_axisym_mw_analysis_real 
<inputFitsFile> <lambda> <J_0> 
<multiresFlag> 

>> . /bin/ s21et_axisym_mw_synthesis_real 

<outputRoot> <lambda> <J_0> <bandLimit> 



Example 1: Performing the forward (analysis) and inverse 
(synthesis) wavelet transform of a real signal (MW sampling) from 
the command line. 



9 Note that the code uses a slightly different notation compared to 
the equations of this article: B refers to the wavelet scaling parameter 



(denoted A herein) and J„ 
Jo herein). 



to the first scale of the transform (denoted 



Example 3: Performing the wavelet transform of a real signal 
(HEALPix sampling) using the Matlab interface. 



4.2. Wavelet transform in Matlab and IDL 

Examples [3] and |4] read real signals on the sphere from FITS 
files, calculate the wavelet coefficients and plot them using a 
Mollweide projection. The first case is a Matlab example 
where the input map is a simulation of the cosmic microwave 
background in the HEALPix sampling. The second case is a 
IDL example where the input map is a topography map of the 
Earth in MW sampling. S2LET ships with versions of these two 
examples in C, Matlab and IDL. 



4.3. Wavelet denoising in c 

Example [5] illustrates the use of the wavelet transform to de- 
noise a signal on the sphere. The input noisy map is a band- 
limited topography map of the Earth in MW sampling at reso- 
lution L = 128. It is read from a FITS file, decomposed into 
wavelet coefficients (for given parameters A and Jq) which are 
then denoised by thresholding. The denoised signal is recon- 
structed from the denoised wavelet coefficients and written to 
a FITS file. 

In this example we consider a noisy signal y = s+n e L 2 (S 2 ), 
where the signal of interest s e L 2 (S 2 ) is contaminated with 
noise n =s L 2 (S 2 ). We consider zero-mean white Gaussian noise 
on the sphere, where the variance of the harmonic coefficients of 
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; Example: Wavelet transform in IDL 
lambda = 3 
JO = 2 

; Read a real MW map from a FITS file 
file = ' data/earth_tomo_mw_12 8 . f its ' 
f = s21et_mw_read_real_map ( f ile ) 
L = s21et_get_mw_bandlimit ( f ) 
Jmax = s21et_j_max (L, lambda) 

; Perform the wavelet transform 

f_wav = s21et_axisym__mw_wav_analysis_real 

(f, lambda, JO) 
f_rec = s21et_axisym_mw_wav_synthesis_real 

( f_wav) 

; Plot the map and the wavelet coefficients 
ns = ceil ( sqrt ( 3+ Jmax- JO ) ) 
!P .MULTI= [ 0, ns, ns] 
s21et_mw_plot_mollweide, f_rec, 

title=' Band-limited map' 
s21et_mw__plot_mollweide, f_wav. seal, 

title=' Scaling map' 
for j=0, Jmax-JO do begin 

s21et_mw__plot_mollweide, f_wav. (j), 

title=' Wavelet map ' tstrtrim ( j+1 , 2 ) +' 
on ' +st rtrim ( Jmax- J0+1 , 2 ) 

endf or 

!P .MULTI=0 



Example 4: Performing the wavelet transform of a real signal (MW 
sampling) using the IDL interface. 



the noise is specified by 

E(|« fm | 2 ) = (r\ W?,m. (17) 

A simple way to evaluate the fidelity of the observed signal y is 
through the signal-to-noise ratio (SNR), which we define on the 
sphere by 



SNR(y) S 101og 10 - ^, 

\\y- s\\i 



where the signal energy is defined by 

\\y\\l = <y|y> = jf d£l(aj)\y(ajf 



(18) 



(19) 



I'm 



We seek a denoised version of y, denoted by d e L 2 (S 2 ), with 
large SNR(c/) so that d isolates the informative signal s. When 
taking the wavelet transform of the noisy signal y, one expects 
the energy of the informative part to be concentrated in a small 
number of wavelet coefficients while the noise energy should be 
spread over various wavelet scales. Since the transform is linear, 
the wavelet coefficients of the j-th scale are simply given by the 
sum of the individual contributions: 



Y j (a)) = S J \a))+N J (oj), 



(20) 



where capital letters denote the wavelet coefficients, i.e. 
Y> sy-k y V j ,S i = s * *F and N j s n ★ *F. For the zero-mean 
white Gaussian noise defined by Eqn. [T7] the noise in wavelet 
space is also zero-mean and Gaussian, with variance 

E(\N J (oj)\ 2 ) = CT 2 Y J \y J m \ 2 ={cr j f. 




(a) Band-limited signal 




(b) Noisy signal with SNR(y) = 11.8dB 




(c) Denoised signal with SNR(<i) = 14.6dB 
Fig. 4: Scale-discretised wavelet denoising by hard-thresholding. 

Denoising is performed by hard-thresholding the wavelet coef- 
ficients F 7 , where the threshold is taken as = 3cr ] . The de- 
noised wavelet coefficients D J = d * are thus given by 



D j (aj) 



0, if YHu>) < THoj) 
Y-i(a>), otherwise 



(21) 



The denoised signal d e L 2 (S 2 ) is reconstructed from its wavelet 
coefficients D' and the scaling coefficients of y, which are not 
thresholded. The denoising procedure outlined above is partic- 
ularly simple and more sophisticated denoising strategies can 
be developed; we adopt this simple denoising strategy merely 
to illustrate the use of the S2LET code. In this example we 
perform the wavelet transform with parameters A — 2 and 
J = 0. For a noisy signal y with SNR(y) = 11.8dB, the 
scale-discretised wavelet denoising recovers a denoised signal 
d with SNR(d) = 14.6dB. The initial, noisy and denoised maps 
are shown in Figure [4] 

5. Summary 

In the era of precision astrophysics and cosmology, large and 
complex data-sets on the sphere must be analysed at high preci- 
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// Example: Wavelet denoising in C 
int lambda = 2, JO = 0; 

// Read a real MW map from a FITS file 
char inputfile [100] = 
double *f; 

int L = s21et_f it s_read_mw_bandlimit ( f ile ) ; 
s21et_mw_allocate_real ( &f , L) ; 
s21et_f its_read_mw_map ( f , file, L) ; 

// Perform multiresolution wavelet analysis 
double *f_wav, *f_scal; 

s21et_axisym_mw_allocate_f_wav_mult ires_real 
(&f_wav, &f_scal, lambda, L, JO); 

s21et_axisym_mw_wav_analysis_multires_real 
(f_wav, f_scal, g, lambda, L, JO ) ; 

// Threshold the wavelets with a noise model 
s21et_axisym_wav_hardthreshold_multires_real 
(f_wav, threshold, lambda, L, JO ) ; 

// Reconstruct the denoised signal 
double * f_denoised; 

s21et_mw_allocate_real (&f_denoised, L) ; 
s21et_axisym_mw_wav_synthesis_multires_real 
(f_denoised, f_wav, f_scal, lambda, L, JO) ; 

// Write the denoised signal 
char outputf ile [100] = 

s21et_f its_write_mw_map (outf ile, f_denoised, L) ; 



Example 5: Denoising a real signal (MW sampling) in C through 
hard-thresholding of the wavelet coefficients. 

sion in order to confront accurate theoretical predictions. Scale- 
discretised wavelets are a powerful analysis technique where 
spatially localised, scale-dependent signal features of interest 
can be extracted and analysed. Combined with a sampling the- 
orem, this framework leads to an exact multiresolution wavelet 
analysis, where signals on the sphere can be reconstructed from 
their scaling and wavelet coefficients exactly. 

We have described S2LET, a fast and robust implementa- 
tion of the scale-discretised wavelet transform. Although the first 
public release of S2LET is restricted to axisymmetric wavelets, 
the generalisation to directional, steerable wavelets will be made 
available in a future release. The core numerical routines of 
S2LET are written in C and have interfaces in Mat lab, IDL and 
Java. Both MW and HEALPix pixelisation schemes are sup- 
ported. In this article we have presented a number of examples 
to illustrate the ease of use of S2LET for performing wavelet 
transform of real signals stored as FITS files and to plot scaling 
and wavelet coefficients on Mollweide projections of the sphere. 
We have also detailed a denoising example where denoising is 
performed through simple hard-thresholding in wavelet space. 
Although only a simple denoising strategy was performed to il- 
lustrate the use of the S2LET code, it nevertheless performed 
very well, highlighting the effectiveness of the scale-discretised 
wavelet transform on the sphere. 
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